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Direct numerical simulation of 
incompressible axisymmetric flows 


By Patrick Loulou 1 


1. Motivation and objectives 

In the present work, we propose to conduct direct numerical simulations (DNS) 
of incompressible turbulent axisymmetric jets and wakes. The objectives of the 
study are to understand the fundamental behavior of axisymmetric jets and wakes, 
which are perhaps the most technologically relevant free shear flows (e.g. combuster 
injectors, propulsion jet). Among the data to be generated are various statistical 
quantities of importance in turbulence modeling, like the mean velocity, turbulent 
stresses, and all the terms in the Reynolds-stress balance equations. In addition, 
we will be interested in the evolution of large-scale structures that are common in 
free shear flows. 

The axisymmetric jet or wake is also a good problem in which to try the newly 
developed b-spline numerical method. Using b-splines as interpolating functions in 
the non-periodic direction offers many advantages. B-splines have local support, 
which leads to sparse matrices that can be efficiently stored and solved. Also, they 
offer spectral-like accuracy and are continuous, where 0 is the order of the 

spline used; this means that derivatives of the velocity such as the vorticity are 
smoothly and accurately represented. For purposes of validation against existing 
results, the present code will also be able to simulate internal flows (ones that require 
a no-slip boundary condition). Implementation of no-slip boundary condition is 
trivial in the context of the b-splines. 

2. Accomplishments 

To simulate these flows, we follow the procedure described in Moser et al (1983) 
and Leonard et al (1982), with b-splines replacing the Jacobi or Tchebyshev poly- 
nomials used in the radial direction. 

2.1 Navier- Stokes equation 

The starting point for the simulations are the incompressible Navier- Stokes equa- 
tions: 

(la) 

V ■ u = 0 (16) 

Let v be our numerical approximation to u, which will consist of a truncated 
expansion in terms of divergence- free vector functions (i.e. V-v = 0). Furthermore, 
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let £ be any vector function representable by another finite set of divergence-free 
vector functions ( V *£ = 0), satisfying f = 0 on dG , the boundary of the domain. By 
substituting v for u in (la) and using the standard weighted residual technique with 
f as the weight functions, we obtain the discrete weak form of the Navier-Stokes 
equations. 


Ia ( ^ dV ~Ia i ' VXUdV = ~Tjj AviV (2 > 

By enforcing (2) for each £ making up a basis for the weight functions, a coupled 
system of ordinary differential equations for the coefficients in the expansion for v 
are obtained, which can be solved using standard time-advance techniques (see sect. 
3). This formulation has the advantage of automatically satisfying the continuity 
equation (lb) and eliminating the pressure. 

2.2 Velocity representation and vector shape functions 
Given the formulation in (2), ail that remains is to select the basis vectors to 
represent v and £. This is facilitated by approximating the spatially developing jet 
or wake of interest as a time developing flow. In this case, the streamwise (axial or 
z) direction is homogeneous and we can approximate the flow as periodic in z with 
period L z . This and the natural periodicity in 9 allow the following representation 
of v and £: 


v(r, 9,z,t) = EEE (3) 

j k , I 

Sj'm-l>(rAz) = W ( 4 ) 

where 

2 nm 

k z — — — , — J < j < J, —M < m < M 

Due to the continuity constraint, there are only two degrees of freedom associated 
with each Fourier/b-spline mode. It is therefore convenient to divide the expansion 
and weight vectors, u / and wi > , into two distinct classes of vectors (u/+ and u/”) 
and (w/' + and w>”), with coefficients o^ m/ (f) and <*” m/ (t). 

The expansion and weight vectors must be constructed such that they are di- 
vergence free and have the proper regularity properties at r = 0 (Shariff 1993). 
The following vectors meet both these requirements provided that 0 j(r) satisfy the 
appropriate condition at r = 0 (see below): 


«r + | 


i<7/(r)\ _ / ~igi(r) 

U/ + (r;;,fc J ) = ( u e + , ] =Vx ( g,(r) J, U/ _ = V x g t (r) ) (5) 

0 / \ 0 




-*£7i'( r ) 




+ = V* x V* x g r (r) , w r =V’xV'x g r (r ) (6) 


0 
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where V x is the Fourier tranformed curl operator and V* x is its complex conjugate. 
The gi(r) are b-splines expansion functions as described by de Boor (1978). 

The above representation is incomplete when k z = 0; so in this case, the following 
special representation is used: 


u/ 


+ 




w / 


,+ _ 




( 7 ) 

( 8 ) 


And when j = 0, these vectors are incomplete; so for that case ( k z = j = 0), we 
use: 



All these additional vectors are also divergence-free and satisfy the regularity re- 
quirement. 

In order to have the correct regularity property, depending on the azimuthal wave 
number j, gi(r) and some of its derivatives (Shariff 1993) must vanish when r = 0. 
Since this is not automatically satisfied by all the b-splines, some of the coefficients 
(a + and a“) must be zero. In particular, for j > 0 


a~l~ =0, 1 < / < min(0,j) + 1 or Ij is odd, while j4-3</<0 + l 

and 

af = 0, 1 < / < min(0 + 1, | j — 1|) or Ij odd \j — 1| + 2 < / < O + 1 

but with af unconstrained. 

2.3 Boundary conditions 

The present method is designed to treat both no-slip boundaries and potential 
boundaries (free shear flows). 

2.3.1 No-slip boundary condition 

Enforcing the boundary condition that u/(i? 2 ) = 0, where J?2 is the outer edge of 
the domain, requires that 
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gi(r = R 2 ) = 0, g',(r = R 2 ) = 0 


( 11 ) 


But if the gt are the b-splines as defined by de Boor, (11) is not satisfied for / = L 
or / = L - 1 (the two functions closest to the boundaries). Therefore, we impose 


a L= a L = a L - 1 = a L - 1 = 0 
However, for the shape functions in (5), (11) also implies 

dug 


dr 


= 0 


( 12 ) 


(13) 


at the boundary, which is too restrictive. To alleviate this, we augment (5) and (6) 
with the additional vectors (Moser et al., 1983): 


0 

u°t — i = ( ~2A:^L_i(r) 
2jg L -i{r) 


ij'9L-i(r ) 

w £-i = 1 + 0 L-i(r) 

0 


( 14 ) 


With this, du$/dr is unconstrained at the boundary. 
i 2.8.2 Potential boundary condition 

When simulating free shear flows, we follow the approach of Corral et al (1993) 
and Sondergaard et al (1994), where it is assumed that the vorticity in the flow is 
confined to a small region (r < JR W ) which is to be computed. In the outer region 
of the flow (r > i?u>), the vorticity is zero, so the velocity is a potential (u = V0 
and A<f> = 0). 

For each Fourier mode, the potential <j> is given by: 

4>{r\j,k z ) ~ I\j(k z r) 

where Kj(x) is the modified Bessel function of the second kind. At the boundary of 
the computational domain, r = R 2 > the following relations are satisfied since 
u is a potential: 


y 

UQ = —U r 

qR2 


ih z 

u 2 = — -u r , where 

9 


_ , K'jjk'R,) 

9 ‘ Kj(k z R 2 ) 


In addition, for the vorticity to be zero at r = R 2 , u r must satisfy 

-K'j(k z R 2 ) 


du T 

7-^r +Ur = o, 


where 


7 = 


(15) 


(16) 


k z K" } (k z R 2 ) 

Given the representation in (5), to satisfy (15) and (16), the coefficients must satisfy: 


a 


[qR 2 


+ 1 


+ a 


L [qR* 


-1 


= 0 


(17) 
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K-i + a L- 1 ) ~ ( a t + a t) + 
where we made use of the following identities: 

9l{R*) = 1, 9 l-ARi) = -9' l{R2), 


19W\ 


1 

7 


1 

R2 


= 0 


j 2 kz_ 
qR\ q 


(18) 


There are three boundary conditions in (15) and (16) but they are redundant 
with the continuity equation at the boundary. Since continuity is built into our 
expansions, only two conditions (17) and (18) are required. 

For the k t = 0 case, the conditions in (15) and (16) are satisfied when: 

c^=0 (19) 


a L - 1 - a t 


1 + 


Rtg'iiR^) 


= 0 


( 20 ) 


3. Future plans 

The only major parts of the code that remain to be implemented are the time 
integration scheme and the non-linear convective term. To time march the equa- 
tions, we propose to use the method of Spalart et al. (1991) which is a mixed 
implicit/ explicit scheme. The linear viscous term is time-marched implicitly using 
a Crank-Nicholson scheme and the non-linear terms are time-marched explicitly 
using a third order Runge-Kutta scheme. To compute the non-linear term, we 
revert back to physical space, take the cross product, integrate exactly by Gauss 
quadratures (doing the integrals exactly takes care of aliasing), and revert back to 
wave-space. 

There are two types of free shear flow simulations that are of interest. First, a 
fully turbulent jet that can be simulated using a turbulent pipe simulation result 
as the initial condition (obtained from the same code). This is similar to the tur- 
bulent mixing layer and wake simulations of Moser and Rogers (1994). Second, a 
transitional jet or wake simulated using different initial conditions, similar to the 
work of Sondergaard et al. (1994) on wakes. 
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